import numpy as np
from scipy.fftpack import fft
import matplotlib.pyplot as plt
import cv2 as cv


def myfft_function(signal, begin, end):
    N = end - begin  # 采样点数
    sample_freq = 20  # 采样频率 20 Hz, 大于两倍的最高频率
    sample_interval = 1 / sample_freq  # 采样间隔s
    signal_len = N * sample_interval  # 信号长度
    t = np.arange(0, signal_len, sample_interval)

    fft_data = fft(signal)

    # 计算幅值
    # 在python的计算方式中，fft结果的直接取模和真实信号的幅值不一样。
    # 对于非直流量的频率，直接取模幅值会扩大N/2倍， 所以需要除了N乘以2。
    # 对于直流量的频率(0Hz)，直接取模幅值会扩大N倍，所以需要除了N。
    fft_amp0 = np.array(np.abs(fft_data) / N * 2)  # 用于计算双边谱
    fft_amp0[0] = 0.5 * fft_amp0[0]
    N_2 = int(N / 2)
    fft_amp = fft_amp0[0:N_2]  # 单边谱

    # 计算相位
    phase_spectrum = np.angle(fft_data, deg=True)  # deg=True表示返回的是角度值，否则是弧度值
    # 将相位从-180~180转换到0~255
    n180 = np.full(N, 180)  # 创建含有N个数的数组，值都为180
    phase_spectrum = (phase_spectrum + n180) * 0.7
    fft_pha = phase_spectrum[0:N_2]

    # 计算频谱的频率轴
    list1 = np.array(range(0, int(N / 2)))
    freq = sample_freq * list1 / N  # 单边谱的频率轴

    return freq, fft_amp, fft_pha, t


save_path = r"E:\studylife\detectflaws\code\python\image\5_26\p60t6"  # 原视频裁切出来的图片的存放地址

column1 = 350  # 无缺陷点坐标（黄色）
row1 = 400
column2 = 400  # 缺陷点坐标（绿色）
row2 = 400
# 读取的图片帧数
map_begin = 1
map_end = 700
# 信号所处频率
freq_origin = 0.45

# 创建列表存储对应帧数视频图像
lists_gray = []
array_gray = []

# 将视频存入array_gray中
for k in range(map_begin, map_end):  # 读取map_quantity帧图片计算fft
    img = cv.imread(f'{save_path}/{k}.jpg', flags=0)  # flags=0读取为灰度图像
    # 为列表赋值
    lists_gray.append(img)
array_gray = np.array(lists_gray)  # 把lists_gray转换成np.array类型方便后续裁切

# 存储各信号灰度值的向量
array_gray1 = array_gray[:, row1, column1]
array_gray2 = array_gray[:, row2, column2]

# 对信号进行傅里叶分解求取频率轴，幅度值，相位值，信号时间
freq1, fft_amp1, fft_pha1, t = myfft_function(array_gray1, map_begin, map_end)
freq2, fft_amp2, fft_pha2, t = myfft_function(array_gray2, map_begin, map_end)

# 打印各信号对应频率下的幅度值和相位值
amp_result1 = np.interp(freq_origin, freq1, fft_amp1)
pha_result1 = np.interp(freq_origin, freq1, fft_pha1)  # 使用np.interp(freq_origin, freq1, pha1)可以获得指定横坐标的纵坐标值
amp_result2 = np.interp(freq_origin, freq2, fft_amp2)
pha_result2 = np.interp(freq_origin, freq2, fft_pha2)  # 使用np.interp(freq_origin, freq1, pha1)可以获得指定横坐标的纵坐标值
print(f"黄色信号锁相调制频率{freq_origin}下的幅度值：", amp_result1)
print(f"绿色信号锁相调制频率{freq_origin}下的幅度值：", amp_result2)

print(f"黄色信号锁相调制频率{freq_origin}下的相位值：", pha_result1)
print(f"绿色信号锁相调制频率{freq_origin}下的相位值：", pha_result2)
print(f"锁相调制频率{freq_origin}下的相位值差：", abs(pha_result2 - pha_result1))

# 绘制结果
plt.figure()
# 原信号
plt.subplot(221)
plt.plot(t, array_gray1, 'y')
plt.plot(t, array_gray2, 'g')
plt.title(' Original signal')
plt.xlabel('t (s)')
plt.ylabel(' Amplitude ')

# 幅值谱
plt.subplot(222)
plt.plot(freq1, fft_amp1, 'y')
plt.plot(freq2, fft_amp2, 'g')
plt.title(' FuPinTu')
# plt.ylim(0, 250)
plt.xlim(0, 1)
plt.xlabel('frequency  (Hz)')
plt.ylabel(' amp ')

# 相位谱
plt.subplot(223)
plt.plot(freq1, fft_pha1, 'y')
plt.plot(freq2, fft_pha2, 'g')
plt.plot([freq_origin, freq_origin], [-10, 280], color='r', linewidth=1, linestyle="--")  # 在对应频率上画出虚线
plt.plot(freq_origin, np.interp(freq_origin, freq1, fft_pha1),
         'r*')  # 使用np.interp(freq_origin, freq1,  pha1)可以获得指定横坐标的纵坐标值
plt.plot(freq_origin, np.interp(freq_origin, freq2, fft_pha2), 'r*')
plt.title(' XiangPinTu')
# plt.ylim(0, 250)
plt.xlim(0, 1)
plt.xlabel('frequency  (Hz)')
plt.ylabel(' pha ')

# 相位差
plt.subplot(224)
plt.plot(freq1, abs(fft_pha1 - fft_pha2), 'b')
plt.plot([freq_origin, freq_origin], [-280, 280], color='r', linewidth=1, linestyle="--")  # 在对应频率上画出虚线
plt.plot(freq_origin, np.interp(freq_origin, freq1, abs(fft_pha2 - fft_pha1)), 'r*')  # 给交点画星号
plt.title(' XiangWeiCha')
# plt.ylim(0, 250)
plt.xlim(0, 1)
plt.xlabel('frequency  (Hz)')
plt.ylabel(' pha_dif ')

plt.show()
